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The dynamical response of edge waves under the influence of self-gravity is examined in an idealized two-dimensional 
model of a proto-stellar disc, characterized in steady state as a rotating vertically infinite cylinder of fluid with constant 
density except for a single density interface at some radius r$. The fluid in basic state is prescribed to rotate with a Keple- 
rian profile ~ r -3 / 2 modified by some additional azimuthal sheared flow. A linear analysis shows that there are two 
azimuthally propagating edge waves, kin to the familiar Rossby waves and surface gravity waves in terrestrial studies, which 
move opposite to one another with respect to the local basic state rotation rate at the interface. Instability only occurs if the 
radial pressure gradient is opposite to that of the density jump (unstably stratified) where self-gravity acts as a wave stabilizer 
irrespective of the stratification of the system. The propagation properties of the waves are discussed in detail in the language 
of vorticity edge waves. The roles of both Boussinesq and non-Boussinesq effects upon the stability and propagation of 
these waves with and without the inclusion of self-gravity are then quantified. The dynamics involved with self-gravity non- 
Boussinesq effect is shown to be a source of vorticity production where there is a jump in the basic state density In addition, 
self-gravity also alters the dynamics via the radial main pressure gradient, which is a Boussinesq effect. Further applications 
of these mechanical insights are presented in the conclusion including the ways in which multiple density jumps or gaps may 
or may not be stable. 
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1. Introduction 

The Atacama Large Millimeter/submillimeter Array (ALMA) has opened up a new window revealing 
the structure of proto-planet accretion discs (pp-discs hereafter). Most notable of the recently reported 
ALMA discoveries is the appearance of asymmetric “peanut” features in the outer (> 80AU) regions 
of several pp-disc systems including SAO 206462, SR 21, LKHa 330 OPH IRS 48 and HD 142527 


(Casassus et a/. 12013) Fukagawa et fl/.|2013| Isella et n/.|2013 

van der Marel et al: 2013 

Perez et al. 

2014 

). The ALMA website recently published a high resolution image of HL Tauri, a 

debris-disc 


phase T-Tauri system, showing a disc with alternating axisymmetric rings and gaps with unambigu¬ 
ous non-axisymmetric features set atop. While planets have yet to be detected in the HL Tauri system, 
the gap-ring features appear to indicate an ongoing process similar to the dynamics responsible for 
the Kirkwood gaps in the Jupiter-Asteroid-belt system, in which an embedded (as yet unseen) planet 
drives gap creation through resonant gravitational forcing of an otherwise structurally uniform Kep- 
lerian disc (A. Parker, private communication). Furthermore, the furthest parts of these same discs are 
likely to be sufficiently cold enough that the densities reached in their local disc midplane’s are cor¬ 
respondingly high enough such that self-gravitational effects should not be ignored. The dynamical 
stability of sheared dynamical structures like those seen in HL Tauri is the main inspiration for this 
work: we seek to develop a theoretical framework that could be helpful in interpreting the dynamical 
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nature of structures emerging in mainly cold pp-disc systems in which self-gravity of the disc gas 
plays an important dynamical role. 

Some pp-disc systems exhibiting peanut features also show indications of severe dust depletion 
within their inner radial locations. It has been suggested that pressure/density enhancements in the 
outer parts of the disc act as particle traps preventing the inner parts of these discs from having 
their dust replenished by accretion from the outer parts. Interestingly, axisymmetric pressure en¬ 
hancements, provided they meet minimum requirements upon their amplitudes and radial girth, are 
known to undergo a barotropic shear transition known as the Rossby Wave Instability (RWI hereafter, 
( Lovelace et «/.||l 999\ Li et al . ||2000 2001)). This instability leads to the creation of coherent anti- 
cylonic vortex structures with wide azimuthal extent. Recent studies have shown that this instability 
can develop even as early as in the infall stage of the disk (Bae et q/. || 2015| >. From both theoretical 
and numerical/computational considerations it is also understood that anticylonic vortex structures in 
pp-discs readily entrain dust through well-known drag-momentum exchange mechanisms ( jBarge and 


Sommeria fl995 1. While the origins of the density enhancement and the causal dilemma posed by the 


above scenario remains unsettled (see recent discussion in Flock et al (2014)), the fact remains that 
rings/gaps and (generally speaking) nearly axisymmetric annular structures are the basic states upon 
which dynamical evolution proceeds in pp-discs. A detailed examination of the dynamical state of 
such structures from a physically mechanistic point of view is therefore of interest to theoreticians 
and observers alike. While the dynamical stability of self-gravitating disc systems is a vast subject 
that we cannot possibly cover here (see the reviews of Armitage (20111), we are interested in ex¬ 
amining the mechanistic nature of stability in self-gravitating zonal flows in accretion discs. Steady 
axisymmetric disc zonal flows arise from the balance between the Coriolis and centrifugal effects and 
the radial dependencies in the pressure and self gravity. Indeed, inspection of the F1L Tauri system 
suggests that the flow within the rings and gaps should exhibit pronounced zonal flow structure since 
the distribution of matter appears nearly axisymmetric, yet radially non-uniform. 

Barotropic zonal flowscan become unstable and the RWI is one example. Barotropic shear insta¬ 
bilities, i.e., instabilities in the absence of baroclinic torque but in the presence of significance rota¬ 
tion may be interpreted within the framework of potential vorticity dynamics and the action of counter 
propagating Rossby waves (Bretherton 1966, Hoskin~ef a/.|198~5] Baines and Mitsuder a| 1 994|[Heifetz 


et a(.||1999[ to n ame just a few), and the RWI has been shown to be consistently interpreted within 
that framework (Umurhan 2010). Potential vorticity [^j (PV) in this context is equal to q/p in which 
q = (V x u) • z is the vorticity in the disc vertical direction and p is the density. The instantaneous 
basic state rotation vector when the disc is in exact rotational balance with no contribution from radial 
pressure support is given by £2fc(r)z with Qk{r) = TT/ f (/'o) (/‘o/r) 3/ 2 where r is radius in a cylindrical 
coordinate system and where ro is a fiducial radius with corresponding rotation rate TLfro). As a gen¬ 
eral rule for axisymmetric basic states, for every radial extremum in the PV gradient there can exist an 
azimuthally propagating Rossby wave radially localized around the corresponding radial extremum 
point. In the simplest configuration in which there is some action-at-a-distance and no possibility of 
critical layers (cf. Marcus et al. (2013)), instability can arise from the resonant interaction of two or 
more Rossby waves (Heifetz et al .\\ 999 1 . For example, in RWI unstable discs with a single axisym¬ 
metric pressure bump, there are two oppositely signed extrema in the radial gradient-PV which can 
each support Rossby waves with wave speeds boosted by the local flow of the disc ( Umurhan|2010[ ). 
Resonant interaction and subsequent instability can occur when the Rossby waves move with equal 


and opposite directions in some reference frame. We also note that the previous studies by Mamat- 


sashvili and Rice (2009) and Mamatsashvili and Chagelishvili| (|2007) have explored the PV profiles 


in order to understand the evolution of the dynamics in PP-discs, 


'The term barotropic in this context refers to lack of variation of the basic shear flow in the vertical direction (i.e, perpendicular to the disk 
mid-plane). The RWI operates whether or not the disturbances are barotropic or baroclinic, but this is different from whether or not the 
basic flow is fully barotropic, fully baroclinic or pseudo-baroclinic. 

2 sometimes known in the astrophysical literature as “vortensity” 
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This perspective of interacting Rossby waves and resonant instability has been extended to include 
other physical effects, e.g. the interaction of resonant Rossby-gravity waves (Hamik et al. |2008| Rabi 


novich et al. |2011 ) and Alfven-Rossby waves in sheared MHD flows (Heifetz et al. |2015|). Generally 


when there is some additional physical effect included in shear flows, baroclinicity begins to take on 
an important dynamical role. In the case of the above-mentioned Rossby-gravity waves, pressure and 
density isolines of the perturbed flow do not coincide (i.e. Vp x Vp ^ 0) and, consequently, a certain 
amount of vorticity creation can occur. Places in a zonal disc flow in which there is a stark density 
gradient can support vorticity waves [^] with this baroclinic “vorticity creation” effect and this, in turn, 
can influence the other vorticity waves of the physical system and consequently change the stability 
character of the whole physical arrangement.A physical perspective developed for a self-gravitating 
sheared system should therefore take this baroclinic effect into account. 

Lovelace and Hohlfeld (j2013|) have extended the examination of the RWI to infinitely thin self- 


gravitating discs and they show that, among other things, transition to instability is very subtle and 
that, if the conditions are right, a localized pressure deficit (as opposed to a pressure bump) can also 
lead to instability. It is the aim of this work to begin constructing a potential vorticity view to use 
in understanding unstable transition in sheared self-gravitating configurations^ Additionally^ careful 
numerical experiments of self-gravitating disc systems performed by Lin and Papaloizou (201 la b) ^] 
examined the 3D disc response of a radial section of a disc supporting a gap cleared out by a perturbing 
planetary body. The numerical experiments show that the resulting azimuthally averaged radial PV 
of their discs exhibit complicated gradient-PV structure, especially in the near vicinity of the two gap 
walls. Interestingly, only one of the gap walls is dynamically stable while the other shows instability 
followed by restricted vortex roll up (see Figure 6 of Lin and Papaloizou (201 lb I). The effort of this 
first of a series of studies is (also) motivated by making mechanistic and rational sense of these results. 

The inclusion of self-gravity in any fluid examination is subtle because it introduces a certain 
amount of ellipticity to the underlying PDE’s. In practice this means having to match the solutions 
developed within the interior of the fluid onto external solutions with appropriate far-held boundary 
conditions. To circumvent the difficulties inherent to this process and simplify the posed mathematical 
problem, most examinations of accretion disc dynamics in the astrophysical literature treat the disc 
as either (i) vertically thin infinitesimally speaking (“razor thin”) with a non-zero vertically integrated 
surface density or (ii) as a vertically infinite cylinder. The former is relatively tractable analytically 
as far held solutions are straightforward to construct while the latter considerably simplihes the sys¬ 
tem because boundary conditions are needed to be specihed in the radial direction only and these, in 
turn, may be represented as parameters of the basic state. Of course, both have severe shortcomings 
if realism is the goal. However, as far as the action of localized edgewaves under the influence of 
self-gravity is concerned, using one or the other of the two approximations ought not contaminate the 
physical content of any result developed forthwith. The model analyzed in this study will be in the 
vertically inhnite cylinder setting, which is basically a two-dimensional system and the simplest of 
the two analytically tractable possibilities. 

In section 2 we formulate the equations of motion in this two-dimensional model setting and de¬ 
velop an evolution equation for the disc vorticity. We further develop both basic states and perturba¬ 
tions. Section 3 concentrates on dynamics where the basic state shows a radial jump in the vorticity 
and density fields. The ensuing interfacial vorticity wave dynamics are deconstructed into the dynam¬ 
ics as arising from standard Rossby edge wave dynamics and the concomitant presence of interfacial 
gravity waves (sections 3.1 and 3.2). The latter of these is further analyzed in terms of Boussinesq and 
non-Boussinesq contributions (sections 3.2.1 and 3.2.2). Most importantly, it is in the latter of these 
subsections where we identify the dynamics arising from self-gravity as being in the Boussinesq and 


3 Although vorticity waves and Rossby waves are often used interchangeably, throughout this manuscript we distinguish between them. We 
consider Rossby waves as associated with vorticity production by the advection of the mean shear, and vorticity waves as associated with 
any vorticity creation. 

4 We note that |Lin and Papaloizou|{201 la[b| consider 2D, razor-thin disks. In the next paragraph we shortly discuss both the razor-thin and 
the infinite cylinders disk approximations. 
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non-Boussinesq category of responses. These deconstructed analysis elements arc reunited in section 
3.3 and shown how they effect the overall stability of the system. It is in this section we show how, 
indeed, self-gravity always has a stabilizing influence irrespective of the radial stratification of the 
basic state. In section 4 we discuss the results and their applications for further analysis. 


2. Formulation 


2.1. Disc equations 

We consider a 2D gaseous disc model whose radial (r) and azimuthal (0) momentum equation com¬ 
ponents are written as, 


Du v 2 l dp dp G grav 
Dt r p dr dr r 2 


M 


(1) 


Dv uv 1 dp dip 

Dt + r pr d6 ^ rdO 


( 2 ) 


The dynamics is described in the inertial frame of the disc where u and v are the radial and azimuthal 
velocities of the flow. We have included in equation (|T]> a model Keplerian potential representing the 
gravitational force of an object with mass M, located at r = 0. Our strategy is to remove this explicit 
potential and rewrite the equations as a deviation around this basic rotationally supported state. Thus 
the basic balance is between the full centripetal acceleration term v 2 /r and the star’s radial force. Thus 
we write v = \’k = rO*, where v* is the basic Keplerian flow, so that 


Jl 

r 


ITrn V T/;:: 


Q.k = 


G grav T/; :: 


( 3 ) 


We then rewrite the azimuthal velocity as this Keplerian base state plus a deviation flow, v i->- v* + v 
where v is henceforth understood to be the deviation azimuthal velocity from this basic Keplerian 
state. As such, equations m are rewritten into the more convenient form as 


Du 

Bt 


-b 2Q. k v I — 


1 dp dip 
p dr dr ’ 


( 4 ) 


Dv (v 1 dr 2 Q.i 


——b u I —I-_. 

Dt \r r dr 

where the Lagrangian time derivative is given as 


1 dp dtp 


prdO rdO' 


dt 


dt dr 


( 5 ) 


D (d , \ d d (v \ d 

Dt = (* + ( u+ ^ r0 )- v ) = d~t +u d;- + yr +ak >^- 


do' 


where u = ur + vO. In this way, the rotationally supported state is built directly into these equations by 
having explicitly removed the potential arising from the central gravitating object. |^] When we refer to 


the rotation rate at some position ro we write £lk(r o) = 


_ / vgrav 




. The disc surface gas density is p, 


p represents the gas pressure, and p is the gravitational potential resulting from the gaseous material 
within the disc, satisfying the Poisson’s relation: 


V 2 p = 4 JtGgrayp. 


( 6 ) 


5 In an exact treatment of a three-dimensional system rotating around a central potential, it is a fundamental fact (e.g. |Kippenhahn et al.\ 
< | 1990^ ) that if the equation of state is assumed to be barotropic then it follows that the basic rotational state will necessarily be indepen- 
dent of the vertical coordinate. The model potential considered in equations CD while not exactly correct as far as a real disc is concerned, 
is meant to be a facsimile of the aforementioned fact. 
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We assume further that the gas is incompressible^ hence 

1 d 1 dv 

v . u= __ M+ __ = o , 

and that, consequently, the local density of the gas is an advected passive tracer: 


( 7 ) 


Dp 

Dt 


= 0 . 


(8) 


Equations (j4|8]) form a complete set to describe the evolution of the disc variables (u,v,p,p,(j>). In 
this paper however we focus on interfacial vorticity waves, thus we transform Q and Q into the 
vorticity equation: 


§- t [(V x u)-z] +u d -^- = ^(Vp x Vp) z 


( 9 ) 


or more explicitly after moving the Coriolis-like term onto the righthand side of the above expression, 


Dq dQ k{ 

- = —U - 

Dt dr 


;e P 


1 

rp 2 


dp dp 
dO dr 


dp dp 
dr dO 


( 10 ) 


where q~ ] -f (rv) — 1is the relative vorticity perpendicular to the disc plane. In arriving to equa¬ 
tion (T0| ) we have made explicit use of equation ([7]). it should be noted that as the flow is incompress¬ 
ible, there is essentially no difference between the vorticity and the vortensity equation. We observe 
that generation of vorticity q, arises from radial advection of the mean disc vorticity <2kep associated 
with the Keplerian shear (Fig. 1 1 (a)), 


^ _1 d(r 2 Q. k ) 

t^kep — -v 

r dr 


(ii) 


and the baroclinic torque (Vp x Vp) • z, whose two terms generate shear and hence vorticity (Figs. 


1(b) & 1(c) i. When isobars and isopicnals intersect, the same pressure gradient exerts larger accelera¬ 


tion on the gas whose density is smaller (simply because the pressure gradient force is — ^Vp). 


2.2. Basic state and linearized dynamics 


Equation set (|4]{8]) admits an axisymmetric steady-state solutions (denoted by bars), satisfying the 
balance of forces: 


1 dp 
p dr 




( 12 ) 


with zero radial velocity (u = 0). This equation is supplemented by the solution of the Poisson equation 
in steady state 

1 d d(j) 

-^-rf- =47rG grav p. (13) 

r dr dr 

together with appropriate boundary conditions. The two basic state equations are characterized by 
four functions of radius, v, p, p and <$>, which means that the model as presented is unconstrained 
and, as a consequence, any basic state flow profile may be specified freely by two arbitrary functions 
of radius. As we are concerned with the dynamics localized primarily in the region of some radial 


'The inclusion of self-gravity effects both compressible and incompressible disturbances in a flow. ?or gravitationally unstable disks it 
is understood that compressible/acoustic dynamics lead to instability. ?e are therefore careful to consider and apply the results of this 
incompressible analysis only to disk conditions where Q > 1, since compressible modes are not known to be unstable when that is the 
case.?urthermore, if the time scales of the dynamics are long enough (compared to the sound crossing times across the relevant scales of the 
system) we may safely consider incompressible/anelastic modeling and examine the effects of SG in these types of modes in isolation. ?t is 
then true that the effects of SG are communicated via the deformation of the interfaces - and it is this process we seek to better understand. 
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position ro, the values of particular basic state quantities evaluated there, like v and (which turn 
out to be important in what follows), can be and are treated henceforth as parameters of the system. 


Linearization of the vorticity equation ( fToj ) with respect to this basic state reveals: 

P L q' _ ,d_ f 1 d(r 2 P. k ) \ 1 

Dt dr r dr J rp 2 


dp dp dp dp 
dO dr dr d6 


where q = 7 f (rv) is the basic state vorticity, small perturbations from the basic state are denoted 

by primes, and ^ = jy + (7 + Qk ) js is the linearized Lagrangian time derivative. Hence, vorticity 
perturbation may be generated by radial advection of the mean ambient (Keplerian plus basic state) 
vorticity 


Q q Qkep q 


1 d(r 2 Q. k ) 
r dr ’ 


and by the linearized action of the baroclinic torque, where its second term is neglected when as¬ 
suming a Boussinesq flow. The baroclinic torque can be written in terms of the radial displacement 
perturbation (u' = 757 q r ) when noting that linearization of the incompressible continuity equation 
([ 8 ]) implies: 


D L p' _ ,dp 
— w _ 
Dt dr 


udp_ 

^ dr ’ 


(15) 


that is to say that all small density perturbations result from radial advection of the basic state density. 
Substitute ( fT5| ) into ( fT4| ) gives: 


P L q' = u ,dQ 1 dp d 

Dt dr rp 2 dr d6 




(16) 


In barotropic incompressible flow material conservation of density implies material conservation of 
pressure and hence, similar to (115 1 , p'i >arotwp i c = Therefore, the sum of the terms in the squared 


brackets represent the residual oaroclinic pressure perturbation. 


3. Interfacial vorticity wave dynamics 


As we mentioned in the Introduction, both numerical experiments and observations suggest that sharp 
gradients of vorticity and density emerge in different stages of the disc evolution. Such shaip gradients 
(either in density or vorticity) will support the generation of interfacial (edge) waves. Here we wish 
to focus on their basic mechanism of propagation and growth and, hence for simplicity, we assume 
strong gradients existing in both the density and vorticity only and replace the corresponding gradient 
terms appearing by 5-functions, i.e., 


dQ 

ih =m{r - 


ro), ^- = ApS(r-ro). (17) 

In practice we are assuming that there is no vorticity production anywhere except for r = ro, This is a 


valid assumption only if || >> | 


therefore we expected that in these cases the 


dynamics can be reasonably represented by the interfacial waves approximation (cf. analogy to MHD 
shear flow Heifetz et al. ( 2015 ])). 

Equation (14 1 implies then that the generated vorticity perturbations must have 5-function structures 
as well. Thus for a wave-like solution the vorticity perturbation takes the form: 

q\r,0,t)=qo8{r-ro)e^ e -^ 


(18) 
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(where the zero subscript indicates evaluation at ro and the azimuthal wavenumber m is taken to 
be positive). Since the gas is assumed incompressible, we can use ([7]) to introduce the perturbation 
streamfunction \/ satisfying: 


/_ 1 dyf' 

U r dO ’ 


/ 

V = 


at/ 


so that q' = V 2 t/. Thus in order for t/ to satisfy (18), it should be in the form: 

x/ = q 0 G(r,rt ) ,m)^ me - at \ 


(19) 

( 20 ) 


where the Green function G (not to be confused with the Gravitational constant denoted by G gmv ) is: 


G(r,ro,m) 


ro 

2m 


fe) ", ifr>r 0 

(^) m , if r — r ° 


( 21 ) 


satisfying V 2 G = \.f.(r^j) — (y) 2 G = 8(r — ro), together with the boundary conditions G(r = 
0,°°) = 0. Note that while u' is continuous, v' is discontinuous at ro and this infinite shear corre¬ 
sponds to the 5-function of the vorticity perturbation there. An example for the structure of t/ is 
shown in Fig. [2] 

We assume a wavelike solution for all the perturbed values in the form of f = f(r. ro.iri)e^ ine <ot \ 
but before solving the general dispersion relation for equation ( fTflj ) we will split the problem to the 
dynamics arises from the existence of the mean vorticity and the mean density gradients. 


3.1. Interfacial Rossby waves 


Only the basic state shear may account to sharp vorticity gradients. In the absence of mean density 
gradient equation ( fTflj ), with the basic state of ( fI7| ), becomes: 

D L q' 


Dt 


= -u A QS(r-ro), 


which together with (p~8||2Tj), provides the dispersion relation: 




m 


-nj = 


A Q 

2m 




ro 


( 22 ) 


(23) 


where Q. 0 is the local rotation rate of the disc at position ro. Equation (22 1 also implies that qo = 


—AQ^o, that is to say, that all variations in the vorticity result from advection of the mean vorticity 
(Fig. 1 1 (a)| ) at the interface of ro. The result is an interfacial Rossby wave which propagates (relative 
to the mean flow) in a direction that keeps the lower mean vorticity to its left. Hence if, for instance, 
the mean flow and the mean vorticity gradient are of opposite signs the waves propagate counter the 
mean flow. The mechanism of propagation is such that vorticity anomalies induce circulation that both 
undulates the interface and advects the mean vorticity to generate fresh vorticity anomalies in concert 
with the undulated interface (Fig. 3(a) & 3(b)] ). These Rossby waves are the interfacial versions of the 
plane-like Rossby waves (Rossby et ai 1939) that are generic in the mid-latitudinal atmosphere due 
to the meridional gradient of the Coriolis force there. 


The dispersion relation of (23) indicates that the waves are neutral on their own. However two 


remote interfacial Rossby waves can generate instability (e.g. Baines and Mitsudera ( j!994j ), [Heifetz 
et al. | ([1999])) by “action-at-a-distance” resonance if the necessary conditions for shear instability are 


satisfied, i.e., that each wave propagates counter the mean flow (Fjprtoft 19531 and that the radial 
mean vorticity gradients at the two interfaces are of opposite sign ( Rayleigh||1880 l. The essence of 
this resonance mechanism is illustrated in (Fig. |4(b)[) and has been applied in order to explain the 


RWI in Umurhan (2010). When self-gravity is included, instability will manifest itself along similar 
lines with the same main requirement - that there are at least two separate interfaces able to resonate 
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with one another. Self-gravity as a dynamical influence only becomes relevant when there is at least 
one density gradient placed somewhere in the flow. The dynamical mechanisms resulting from such a 
single density gradient remains to be elucidated and this is what we do next. 


3.2. Interfacial gravity waves 

In the absence of mean vorticity gradient together with the basic state expressed in equation ( [T7] ), 
equation (|T6]) becomes: 


D L q' d 


Dt 


de 


\,dp , 

S ^~+P 

dr 


^E§(r-r 0 ), 

rp z 


(24) 


where po is evaluated as the averaged mean density across the jump . We first obtain the solution for an 
interfacial gravity wave under the Boussinesq approximation which neglects the pressure perturbation 
in the squared brackets. 

3.2.1. Boussinesq waves 

f (j uu. ui^jjiuv.v.nn.111 ULjuauuii i 1 “ 

equation there yields: 


At the position r = ro the cross-stream displacement equation it = together with the vorticity 


q 0 = -2m (— -Q 0 ) fo, 

(- 

~&o 

70 = - 

\m / 




and the corresponding dispersion relation is: 




Bous o 

i 

Ap / 


=t= 

m 

2 rm 

P v 

kP dr ). 


Ap dp 
rp 


( -*2 d r ) 


ro 


1/2 


ro 


(25) 


(26) 


The wave is neutral when the signs of Ap and ^ agree. For instance, if the mean flow is relatively 

weak, (12) is close to hydrostatic balance and under stable stratification (Ap < 0) ^ < 0. The wave 

propagation mechanism for such scenario is illustrated in Figs. 3(c,d). For clockwise propagation 

'OF 


relative to the mean flow, 


Bous 

m 


— D. 0 I <0, equation set (251 determines that q' and t,' are in phase, 


hence the radial velocity field (indicated by the radial arrows) induced by the vorticity perturbation 
shifts the displacement in the clockwise direction. Since Ap < 0, ( [15] ) indicates that positive radial 
displacements corresponds to positive density anomalies, thus for this wave q' and p' are in phase 
as well. Now the Boussinesq part of the baroclinic torque results from the first term at the squared 
brackets of ( fj~4| ) and according to the mechanism described in Fig. 2(a) this will shift the vorticity 
anomalies clockwise, in concert with the radial displacement. For the counter-clockwise solution the 
vorticity and the displacement anomalies are in anti-phase where the displacement and the density 
anomalies remain in phase. 

The unstable stratified case (Ap > 0 & < 0) generates instability of the Rayleigh-Taylor type. 

The wave is advected by the mean flow vq and grows with a rate given by 


m 
2r 



(27) 


The perturbation potential vorticity q' is shifted by a quarter of a wavelength counter-clockwise to q' 
and p' and q' are in anti-phase (Fig. [4(a) l. The resonance amplification results from the fact that in 
this configuration u' is in phase with q', and the Boussinesq baroclinic torque is in phase with q'. 

We thus see that both the propagation rate and the instability mechanism of Boussinesq gravity 
waves depend heavily on the sign of the mean radial pressure gradient force. Flence, it is crucial to 
include the contribution of the mean self gravity in ( fT2j ) to the calculations, if the mean potential 
gravity gradient is comparable in magnitude to the Coriolis and the centrifugal forces. 
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3.2.2. Non-Boussinesq waves 


Considering only the non-Boussinesq component of the baroclinic torque, equation (24) yields: 
D L q' dp' Ap 


d, = S0 7p S(r - ro ’ : 




\m 


rp- 


Po ■ 


(28) 


ro 


In order to evaluate po we first linearize the azimuthal component of the momentum equation ([5]) to 
obtain: 


D.v' - . 1 d (p' . 

Dt ^ r d6 \p Y 


(29) 


and note that since v' is discontinuous and changes sign at ro, vo must remains zero at all times. Hence 
at ro the wave satisfies the angular momentum balance between the radial flux of the absolute vorticity 
and the perturbation azimuthal torque: 


r »2» 4 = 


ro 


which with the aid of (fl9|) becomes: 


Po = Po(QoV4) + Po)- 


(30) 


(31) 


In the absence of self gravity perturbations (po = 0), we can substitute ipo directly in (24) to obtain, 
together with the displacement stream function relation, the dispersion relation: 


7i~ a " + { 2 ^ a 


( Ap 


(32) 


ro 


Consider for instance the stably stratified case where Ap < 0 and Q > 0 (the sign of the absolute 
vorticity matches the sign on the Keplerian one) so that the wave is neutral and propagates in the 
clockwise direction relative to the mean flow (Fig. |3(e)| ). Since the stream-function and the vortic¬ 
ity perturbations are in anti phase at ro, equation ( f3T| ) indicates that the pressure and the vorticity 
perturbations are in anti-phase there. Then the non-Boussinesq component of the baroclinic torque 
(proportional to — ( 57 ^) will generates fresh vorticity anomalies quarter of wavelength in the clock¬ 
wise direction. Since the displacement Q should be translated in concert with the vorticity this implies 
that it should be in phase with the vorticity anomaly. We note that this w ave c annot become unstable, 
even in unstably stratified setup. For counterclockwise propagation (Fig. 3(f)) Ap < 0 and Q < 0 is 
in anti phase with q' and //. 

In order to isolate the contribution of the self gravity perturbation we may assume a quasi- 
hydorstatic balance; po = popo, in equation |3T| ). Combining then ^ with (18) and the Green’s 
function solution gives V 2 p' = —4nG grav ApS(r — ro)^', hence 


i(m0—(Ot) 


P' = -4nG grav ApoG{r,ro,m)loe 

^ Po = POpO = (^—^ r oGg r avPoApo ) <§0- 


(33) 


The last relation manifests the quasi-hydrostatic balance; if, for instance Apo < 0, an outward radial 
displacement shifts mass outward and therefore decreases both the gravitational potential and the 
pressure perturbation that balances the gravity fluctuation. Substitute ( |33| ) in ( |24| ), together with the 
generic relation obtained from the cross-stream displacement equation u' = we obtain: 


qo = -2m(—-£l 0 )%o, 

\m 


[CO 


2n 


1 ) 9 ° — KJ grav 

V m / m 


(Ap ) 21 

p 


Zo, 


(34) 


ro 
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with the corresponding dispersion relation: 


CO _ 

— =ft 
m 


'o ■ 


—r 
m A 


(Ap) 2 


21 1/2 


grav 


( 35 ) 


'll 


It is not surprising that quasi-hydrostatic fluctuations cannot generate instability. Consider again a 
relative clockwise propagation (Fig. 3(g) l, then when Ap < 0 at ro, (c/■ ,p') are in phase with each 

other and in anti phase with (p r ,p') there, ensuring the coherent propagation of the wave. Even in the 
unstable stratified case (Ap 0 > 0), the wave remains neutral and propagates clockwise if (c/. c,'. p' ,p') 
are in phase with each other and in anti phase with p' at ro. For counterclockwise propagation (Fig. 


3(h) i is in anti phase with (</, p',p') 


3.3. Mixed Rossby-gravity interfacial waves 


We can now combine the effects of all four components (Rossby, Boussinesq, the non-Boussineq 
vorticity flux and the non-Boussineq self gravity perturbation) to derive the general interfacial dis¬ 
persion relation of ( p~6| ). for the sharp gradients of equation < [T7| ). We apply the same strategy applied 
previously following all the substitutions we made in subsections (3.1,2), followed by simplifying 
equations ( fT9j[2T ) to derive an expression for uq, and further combining and simplifying equations 
(30, 32) to obtain an expression of po- These two terms are both expressed in terms of the mean flow 
properties and (qo, /jo). This process results the generalized dispersion relation: 


--a„) = T 

/ 4/77 


m 


AQ + ^-Q 


ro 




l 2 


Ap 

P 


2 mr 


1 dp 
p dr 


^ + 


KGgravAp 


ro 


J grav•- 
m- 


1/2 


ro 


together with the generic vorticity displacement ratio: 


co 


q Q = -2/7?-£2 () § 


m 


(36) 


(37) 


All the restoring mechanisms which, on their own, support neutral propagating waves now acting 
together to stabilize the only Rayleigh-Taylor like Boussinesq instability which is obtained when the 
product Ap • fr < Therefore, the necessary criterion for instability is the same as for standard 
Rayleigh-Taylor instability 


» 3 ir 


< 0 . 


ro 


With the inclusion of self-gravity, the sufficient condition for instability for a single edge wave is 

2 r 0 




< -- 


ro 


m 


( pAQ + ApQ\ 2 


V 


4 




4" ftGg rav pAp 


(38) 


r o 


The total composite expression within the square brackets on the right hand side of equation (381 is 
always greater than zero, thus the righthand side is always less than zero implying that these non- 
Boussinesq terms act in concert to stabilize the familiar Rayleigh-Taylor instability. 


4. Discussion 

This single-interface setup and analysis has been meant to be one that showcases the propagation and 
instability of individual vorticity waves under the influence of mean vorticity and density gradients in 
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the presence of self-gravity. The direction of propagation and the ability of the wave to be unstable 
boils down to the phasing relationship between the vorticity and the radial displacement anomalies. 
When the vorticity and the displacement anomalies are in phase (antiphase), the result will be a clock¬ 
wise (anti-clockwise) propagating wave as in fig 2 (a)|2(b)| , when the displacement is in a quarter of a 
wavelength behind the vorticity, the wave will become unstable as in fig |2(c) It is not surprising that 
where there was only one edge wave present in the absence of gravitational dynamics, upon its inclu¬ 
sion (both through the global pressure gradient and self-gravity) an edge can support two so-called 
kemal-gravity waves ( [Harnik et al. |2008 ). It is also not surprising that one recovers the necessary 
condition for Rayleigh-Taylor instability if the flow is unstably stratified in the classical sense, i.e. if 

AP^<0. 


However, it is surprising that the stability criterion found in equation ( [38] ) implies that in the single 
interface model kernel gravity waves are directly stabilized by self-gravity and its stabilizing role 
is unaffected by the stratification of the background state. One could easily imagine there being a 
location in a disc where the local background pressure gradient is identically zero (because of, say the 
global radial distribution of the gas material). In that case the effective radial component of gravity is 
zero and there ought not be any Rayleigh-Taylor like dynamics present. Under such circumstances, if 
the density jump across the interface is non-zero, the consequent activation of self-gravity works to 
stabilize disturbances and its dynamical effect is to alter the vorticity wave propagation speed only. 

It is also worth noting that the role of each of the following effects, self-gravity, Rossby and the 
non-Bousinessq arc to stabilize the flow. This influence is diminished with decreasing wavelength 
distributions leaving short wavelength unaffected, as can be seen from equation ([36]). The result ap¬ 
pears consistent with intuition as the dynamical influence of self-gravity weakens at small scales. 

Interestingly, self gravity can also indirectly influence the onset of the Rayleigh-Taylor instability 
through the basic state (Eq. 12) and possibly render the flow to be unstable. For example, in the case 
that the self gravity term dominates over the inertial terms in the RHS of equation ( p~2| ). one can 
conclude that these dynamics, which depends strongly on the basic state radial pressure gradient, is 


depe 

€ 


mainly due to self gravity contribution! 

The relevancy of self gravity to the dynamics is usually determined by the value of the Toomre 
parameter which is defined as Qj = K Q kC \ for Keplerian disk; where C s is the local sound speed and 
Z is the vertically integrated density. This parameter was originally derived for razor thin disk and 
it was showed that the disk is stable against gravitational instability for axisymmetric perturbation 
if Qt > 1- It can be shown that this parameter is related to 


h M* 


(Armitage 


2010 ) (where h is the 


1 Mrfi s k _ ^ 

pressure scale height), which is a more intuitive interpretation of the gravitational instability. It states 
that the favorable conditions for the onset of the instability are at places in the disk which are flat 
(cold) and heavy. Figure [3] shows the Toomre parameter as a function of the natural logarithm of the 
radius. It is easily seen that the expected value of Qt is highly uncertain and depends on the chosen 


model (which determines -) and the ratio of 


M t 

Mdis 


These quantities are unfortunately unknown both 


theoretically and observationally. 

Moreover, many studies have shown that disks can be gravitationally unstable for values of Qj 
slightly higher then the classical values. In addition, recent studies have shown that self gravity is 
important to the overall dynamics even if the value of Qt is significantly higher then unity. It is known 


that self gravity has a major effect on the RWI instability in many stages of the dynamics (Fin and 


Pap aloizou|201 lb Fovelace and Hohlfeld|2013[ ) , and in spiral arms generation (Fin and Papaloizou 
2011a Goldreich and Tremaine||1979). A form of the Q-parameter in this coordinate system already 


appears in equation (36) and is given by 


jLin and Papaloizou| l[201 lb) also report opposite contributions to SG in the perturbed and basic states in their study. However, the system 
they considered is that of two walls and the interaction could be much more complicated and should be examined in this context. 
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P^o 


XG g rav&P~ 


It is straightforward to calculate a critical value of for the Rayleigh-Taylor instability discussed 
earlier by equating the two terms in the brackets of the RHS in equation (38). It can be seen that £}p 
depends on the relative sizes of the density and vorticity jumps. Nonetheless, one can deduce that 
SG has an important effect on Rayleigh-Taylor instability even for very high £?r values. We therefore 
conclude that for a single interface, instability occurs only if the flow is unstably stratified with respect 
to the sign of the density jump (Fig. 4(a)), and self-gravitational effects should be considered even for 
Toomre parameter values exceeding one. In the absence of instability self-gravity alters the speed of 
propagation of the two kernel gravity waves present along the individual interface. 

In a follow-up study we currently examining the role of both self-gravity and non-Boussinesq ef¬ 
fects in the instability mechanism between well separated edge waves (Fig. |4(bj i, propagating along 
locally stably stratified regions, for different Toomre parameter regimes. 
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FIGURES 





Figure 1. Mechanisms of vorticity generation q (indicated by a counterclockwise dashed circle) in cylindrical geometry, (a) Radial ad- 
vection of the mean disc vorticity (accounted for the Rossby wave mechanism, Section 3.1). The thick arrow represents the radial velocity, 
(b&c) The two components of the baroclinic torque (accounted for the Boussinesq and non-Boussinesq wave mechanisms, respectively. 
Section 3.2,3). Black arrows indicate the direction and relative strength of the velocity induced by the torque. The order of the terms at the 
RHS of equation (7) corresponds to the order of the diagrams. 
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(c) 


Figure 2. Contour plots of the streamfunction i// (16) for an azimuthal wavenumber m = 5, centered at the interface ro = 1. Positive 
values of if/ correspond to negative values of q', hence the red arrows indicate the direction of the radial perturbation velocity = u! 
corresponding to t//. The black wavy line represents the interface displacement at present time, while the red dashed one represents its 
evolution at some time t + At due to the advection of u'. Scenarios in which the vorticity and the displacement are (a) in phase resulting 
clockwise propagation (with respect to the mean flow) (b) in anti-phase resulting counterclockwise propagation (c) in quadrature resulting 
instability. In all figures the Keplerian rotation is counterclockwise. 
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(b) 
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(e) 
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FIGURES 




Figure 3. Vorticity waves propagation mechanisms. The black wavy lines represent present interfacial displacements while the red dashed 
ones represent their evolution (as in Fig. 2). On odd (even) figures the waves are propagating clockwise (counterclockwise), where each 
pair of clockwise and anti clockwise figure corresponds to a different wave propagation mechanism: (a,b) Rossby, (c,d) Boussinesq gravity, 
(e,f) Non-Boussinesq gravity, (g,h) Non-Boussinesq self gravitating waves. In (a)-(h) the background white (gray) color indicate positive 
(negative) signs of cj . The red dashed circular double head arrows indicate the vorticity production. Note that the Keplerian rotation is 
counterclockwise, more details about the different propagation mechanisms appear in Section 3. 
















Geophysical and Astrophysical Fluid Dynamics 


Self Grav' Ros' submitted' refine' 09 ’ 01'16 


FIGURES 


19 



(a) SG modified RayleighTaylor instability 



(b) 2 interfacial Rossby waves action-at-a-distance instability 


Figure 4. Instability mechanism of interfacial vorticity waves, (a) Single interface - when the vorticity (gray circles) and the displacement 
(wavy gray line) are in quadrature radial velocity induced by the vorticity (doubled red arrows) amplify the displacement (dashed red wavy 
line) and the baroclinic torque amplify back the vorticity (red dashed circles) resulting in positive feedback, (b) An illustration of a shear 
flow which can support two counter-propagating vorticity wave resonance. Lines and arrows are as in the previous figures. Each wave 
propagates counter its local mean flow and therefore the waves may become phase locked. Furthermore, the radial velocity induced by 
each wave on the opposed one (indicated by the smaller evanescent red arrows) amplify the other wave’s amplitude. As a result continuous 
mutual amplification of the waves’ amplitudes is obtained. As before, the Keplerian rotation is counterclockwise. 
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FIGURES 



Figure 5. Toomre parameter as a function of the natural logarithm of the radius, for a disk rotating around a \M sun star for radii of 0.4-270 
[AU]. All lines correspond to radiative hydrostatic equilibrium model of |Chiang and Goldreich|(l997| hereinafter CG97), the color indicates 
what was selected for the calculation of for black £ is as in CG97, the same for the green line but with j. adjusted to 0.05 at 1 AU, and 
in the red line this ratio is constant with constant * =0.05. Also the full, dashed and dotted lines corresponds to values of E at 1 AU of 
(1,5,10) * 10 3 [^t] respectively (where 1 * 10 3 [^y] correspond to this disk MMSN). It should be noted that the jumps in side the lines, 
seen especially in the black lines, is due to the fact the CG97 model isn’t continues for J, and is divided into three regions. 













